lib/Mi manual de algoritmos/version_actual/src/geometria/circle_through_3_points.cpp

   1 point center(const point &p, const point &q, const point &r) {
   2     double ax = q.x - p.x;
   3     double ay = q.y - p.y;
   4     double bx = r.x - p.x;
   5     double by = r.y - p.y;
   6     double d = ax*by - bx*ay;
   7
   8     if (cmp(d, 0) == 0) {
   9         printf("Points are collinear!\n");
  10         assert(false);
  11     }
  12
  13     double cx = (q.x + p.x) / 2;
  14     double cy = (q.y + p.y) / 2;
  15     double dx = (r.x + p.x) / 2;
  16     double dy = (r.y + p.y) / 2;
  17
  18     double t1 = bx*dx + by*dy;
  19     double t2 = ax*cx + ay*cy;
  20
  21     double x = (by*t2 - ay*t1) / d;
  22     double y = (ax*t1 - bx*t2) / d;
  23
  24     return point(x, y);
  25 }